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We compare the performance of an UV-filtered Multiboson algorithm, including a global quasi-heatbath update, 
to the standard Hybrid Monte Carlo algorithm in full QCD with two flavours of Wilson fermions. 



1. Introduction 

In this study we consider QCD with two de- 
generate flavours of dynamical Wilson fermions. 
The partition function is given by 

Z = y"nrff^=^,P e-^'""' det(^[[/]+m)2 . (1) 



The fermion determinant which results after inte- 
gration of the fermions leads to non-local interac- 
tions among the gauge links Ux.^ making the cost 
for updating all links grow naively like OiV"^), 
where V is the lattice volume. 

In the Hybrid Monte Carlo (HMC) algorithm, 
which is the standard method of numerical simu- 
lations with dynamical fermions, one uses an aux- 
iliary bosonic pseudo-fermion field (j) to write the 
determinant as 



|det(^ +m)p = / [d(t)''][d(j)] e 



(2) 



and then uses Molecular Dynamics to globally up- 
date the U. 

An alternative approach, which allows the use 
of local algorithms, is the Multiboson (MB) 
method, originally proposed by Liischer A 
polynomial Pn{x) — CnY[k=ii^ ~ ^k) is con- 
structed which approximates 1/x over the whole 
spectrum of {Ip + m). Then PniP + rn) « 
{p + m)~^ and one can write 



I det(^ + m)p « |det"^P„(^ -I- to)^ 



(3) 



In this way the original QCD partition function is 
written approximately in terms of a local action 
so that standard powerful local algorithms (heat- 
bath, overrelaxation) can be used. The system- 
atic error deriving from this approximation can 
be corrected either during the simulation, with 
a global accept-reject step, or by a reweighting 
procedure on physical observables. 

We compare here the efficiency of an improved 
version of the MB algorithm with a "state of the 
art" HMC algorithm, namely the one used by the 
SESAM collaboration [^. The observables used 
for this comparison are the plaquette, represen- 
tative of the smallest (UV) scale, and the topo- 
logical charge, representative of the largest (IR) 
scale features of the gauge field. 
2. Algorithm 

We use the exact, non-hermitian version of the 
MB algorithm |^ with a noisy Metropolis test to 
correct for the polynomial approximation to the 
fermionic determinant. There exist two different, 
complementary strategies for improving the effi- 
ciency of the MB method: 

l.UV-filtering: The first is the choice of an op- 
timal polynomial P„ which keeps n as low as 
possible and still achieves a reasonable Metropo- 
lis acceptance. We follow the procedure of jj] 
inspired by the loop expansion of the fermion 

determinant: det(l — kM) — e X^j^^Tr*'^ x 
det((l-KM) e'^^^"'^']. We group the 
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term e ^j""^^^^ with the gauge action (for 
j < 4 this amounts to a mere shift in /?) and 
find a polynomial approximation to the inverse 

of (1 — kAI) e^^j "^*^ . The parameters have 
to be chosen so that det W ^ 1, where W = 

111(1 - kM ~ Zkl) ■ (1 - kM) ■ e^r=o' . A 
sufficient condition for this is that Wrj ~ i] for 
Gaussian vectors rj. Therefore one takes an equi- 
librium gauge configuration, fixes the coefficients 
aj to some initial value, draws one or more Gaus- 
sian vectors 77 and finds, by quadratic minimiza- 
tion, the roots {z^} which minimize the quantity 
e =11 Wr] — rj |p. The derivatives de/daj are also 
obtained, which allows the optimization of both 
the coefficients {aj} and the zeros {zk\^ following 
ref. |. 

2. Quasi - heatbath: The second is the choice of the 
update algorithm for gauge links and boson fields. 
The coupled dynamics of the system (gauge links 
-f boson fields) is highly non trivial, so that the 
optimal mixture of update algorithms is essen- 
tially determined by numerical experiment |^]. 
For light quark masses one needs a good algo- 
rithm to speed up the IR bosonic modes. To im- 
plement a global heatbath on the bosonic fields 
one must assign 0^^^ = {D — Zk)~^rj where D = 
1 — kM and 77 is a Gaussian random vector. This 
procedure is expensive especially for small quark 
masses. Instead we use a quasi-heatbath consist- 
ing of an approximate inversion plus a Metropolis 
accept-reject step as proposed in Q. i.e. we solve 
approximately {D — Zk)x = [D — Zk)<l)2^^ -\- rj and 
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with r the residual. We combine this global 
update with local overrelaxation of gauge and 
bosonic fields. 

3. Results 

Three representative systems, denoted (A),(B) 
and (C), in Table I, have been studied: medium- 
heavy quarks in a small lattice, light quarks in 
a small lattice, and light quarks in a large lat- 
tice. The HMC simulations used for compar- 
isons of these systems incorporate state-of-the-art 



improvements (but without the recently imple- 
mented SSOR preconditioning). 



Simulation 


A 


B 


C 


Volume 


84 


84 


16^ X 24 


13, n 


5.3,0.158 


5.3,0.165 


5.6,0.156 


'^bosons 


7 


16 


24 


a2 


1.389 


6.066 


4.077 


a4 


4.411 


4.423 


8.789 


A/3 = 192K^a4 


0.527 


0.629 


0.999 


nntia){MB) 


- 3500 


- 64000 


- 27500 


tm{^){HMC) 


- 14000 


- 72000 


- 85000 



Table 1 

Summary of the parameters of our 3 simulations. 
Integrated autocorrelation times are measured in 
units of applications of the Wilson Dirac matrix 
to a vector. 
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Figure 1. Autocorrelation of the plaquette, with 
the UV-filtered MB algorithm (solid line) and the 
HMC algorithm (dashed line) at /3 = 5.3 and n = 
0.158 (simulation A). 

Medium - heavy quarks on a small lattice (Simu- 
lation A): Only seven bosonic fields are required, 
with only one zero Zk devoted to the UV part of 
the Dirac spectrum. As demonstrated in Fig. 1 
the plaquette decorrelates about 4 times faster 
than with HMC. 

Light quarks on a small lattice (Simulation B): 
To see if there is a quark mass below which the 
MB loses its advantage we consider light quarks. 
Several of the IR zeros are now closer to the Dirac 
spectrum boundary. The quasi-global heatbath 
provides an improvement of at least a factor of 5. 
From the plaquette autocorrelation the HMC and 
MB algorithms appear equally efficient. 
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Figure 2. Zeros of the UV-filtered polynomial 
(n = 24, squares), of the non- UV-filtered polyno- 
mial (n — 80, diamonds), and estimated bound- 
ary of the Dirac spectrum (x), (simulation C). 
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Figure 3. Autocorrelation of the plaquette, with 
the UV-filtered MB algorithm (solid line) and the 
SESAM HMC algorithm (dashed line) at /3 = 5.6 
and K — 0.156 (simulation C). 

Light quarks on a large lattice (Simulation C): 
The number of boson fields needed is 24, as com- 
pared with 80 which would have been needed 
with no UV-filtering for a comparable acceptance. 
Fig. 2 shows that 16 of the roots of the UV-filtered 
polynomial are devoted to the IR modes and 8 to 
the UV modes of the Dirac operator. 

The autocorrelation of the plaquette is com- 
pared in Fig. 3 with that obtained by the SESAM 
collaboration using HMC (|. Our MB approach 
is more efficient by a factor '^3. In Fig. 4 we 
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Figure 4. Comparison of topological charge histo- 
ries obtained with HMC and UV-filtered MB al- 
gorithm at /3 = 5.6 and k = 0.156. The common 
scale has been set in terms of equivalent D x v 
mult iplications. 

compare the topological charge histories obtained 
from our simulation and from a sample of SESAM 
configurations j^]. In both cases the same cool- 
ing method was used. Neither simulation is long 
enough to extract a reliable autocorrelation time 
for this observable. Attempts at doing so yield 
roughly equivalent results for both algorithms, as 
the figure already indicates. Therefore, even for 
this global observable, our MB method seems not 
worse than HMC. 

4. Conclusions 

The numerical evidence presented in this study 
shows that the non-hermitian MB algorithm with 
UV-filtering and a global quasi-heatbath of the 
boson fields is a superior alternative to the tra- 
ditional HMC: it decorrelates the plaquette more 
efficiently, and the topological charge equivalently 
well. 
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